Load some packages that we’ll need to use to do these calculations:

library(tidyverse)
library(gifski)
library(ggraph)
library(here)
library(igraph)
library(tnet) # for the closeness function described here: https://toreopsahl.com/2010/03/20/closeness-centrality-in-networks-with-disconnected-components/
source(here("modelFunction_rewiring.R"))

1. Run the model once

# Define parameters
N = 50
socAlpha = 2
socBeta = 4
n.removed = 10
burn.in = 10
recovery = 5
mod00 = -0.4
mod01 = 0.2
mod10 = -0.2
mod11 = 0.4
coefGain = -1 # next step is to make these effects depend on the sociability and current situation of the individual.
coefKeep = -1
modelGraphs <- runModel(N = N, # Nodes in the network
                        socAlpha = socAlpha,
                        socBeta = socBeta,
                        n.removed = n.removed,
                        burn.in = burn.in,
                        recovery = recovery,
                        mod00 = mod00, 
                        mod11 = mod11, 
                        mod10 = mod10, 
                        mod01 = mod01,
                        coefGain = coefGain,
                        coefKeep = coefKeep)$graphs

2. Do individuals have differentiable degrees over time?

df <- lapply(modelGraphs[1:burn.in], degree) %>%
  do.call(rbind, .) %>%
  as.data.frame() %>%
  mutate(timestep = 1:nrow(.)) %>%
  pivot_longer(cols = -timestep, names_to = "id", values_to = "degree") %>%
  group_by(id) %>%
  mutate(initDegree = degree[1]) %>%
  ungroup() %>%
  mutate(id = as.numeric(stringr::str_remove(id, "V")))

df %>%
  ggplot(aes(x = timestep, y = degree, col = initDegree, group = id))+
  geom_line()+
  scale_color_viridis_c()

# yes, we definitely see some individual differences in degree over time. The strength of this effect depends on socAlpha and socBeta.

4. Relationship between removed nodes and network-level measures

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

A few sanity checks:

What about the ratio between the first and second changes? Aka: what percentage of the loss/gain is recovered by the rewiring?

## Warning in log(ratio): NaNs produced

## Warning in log(ratio): NaNs produced

## Warning in log(ratio): NaNs produced
## Warning: Removed 1005 rows containing non-finite values (stat_smooth).
## Warning: Removed 986 rows containing missing values (geom_point).

## Warning in log(ratio): NaNs produced
## Warning in log(ratio): NaNs produced

## Warning in log(ratio): NaNs produced
## Warning: Removed 1005 rows containing non-finite values (stat_smooth).
## Warning: Removed 986 rows containing missing values (geom_point).

## Warning in log(ratio): NaNs produced
## Warning in log(ratio): NaNs produced

## Warning in log(ratio): NaNs produced
## Warning: Removed 1005 rows containing non-finite values (stat_smooth).
## Warning: Removed 986 rows containing missing values (geom_point).

## Warning: Removed 19 rows containing non-finite values (stat_smooth).

## Warning: Removed 19 rows containing non-finite values (stat_smooth).

## Warning: Removed 19 rows containing non-finite values (stat_smooth).